De manera abstracta podemos ver que Los analistas existen para brindar explicaciones y recomendaciones sobre algún asunto de interés. El decisor elije el curso de acción sabiendo que no puede esperar del analista información completa. El gestor se encarga de implementar las decisiones; esa implementación traerá nueva información y el analista vuelve a su trabajo.
Para acabar este curso sintiéndote que vas en camino a ser analista político, te recomiendo seguir estos pasos:
Identificar algún patrón (no aleatorio) en la complejidad del mundo, que afectan los objetivos de las organizaciones: el problema. La identificación del problema depende de cómo entiendes al mundo, por lo que debes haber leído diversas teorías antes de aventurarte a identificar.
Estructurar lo identificado (el problema) en preguntas de investigación.
Revisar la literatura para saber qué respuestas hayan dado otros investigadores al problema planteado.
Plantear explicaciones al problema: hipótesis.
Organizar una estrategia para probar, o invalidar, las hipótesis.
Aplicar la estrategia, entender y validar los resultados. Asegurar la replicabilidad de la estrategia.
Elaborar síntesis de lo actuado; proponer explicaciones a lo encontrado; y enviar recomendaciones.
Hay muchas opciones para la estrategia señalada en el punto 5. La elección de las mismas dependerá de la preparación del analista. En este curso, te daremos conocimientos para que no desestimes el camino basado en datos y/o evidencias empíricas. De ahí que parte de tu proyecto para este curso está basado en desarrollar los puntos anteriores.
Con tu hipótesis (punto 4), debes organizar un archivo de datos (punto 5), donde, al menos, una de las columnas sea la variable que representa el problema observado (la dependiente), por ejemplo: muertos por Covid, resultados de exámenes escolares, etc.; y las otras (las independientes) sean variables sugeridas en la revisión de la literatura como posibles influencias para el comportamiento de la variable dependiente. Refresca estos conceptos con esta lectura.
Con esos datos, prueba tu hipótesis con un análisis de regresión (punto 6). Hay varios pasos a seguir, y lo veremos en las primeras sesiones.
Todos deben colaborar con diversos datos, pues luego debes usarlos para análisis de conglomerados y factorial. Aquí necesitan diversas independientes. Éstas pueden ser dos grupos, uno que represente indicadores similares entre sí (vea el índice de democracia) o disímiles entre sí (vea el índice de desarrollo humano). Esto debe permitirte un mejor análisis de regresión.
Las investigaciones, en general, y las políticas, en particular, deben ser auditables. De ahí que debes organizar todo el trabajo anterior en un repositorio en la “nube”. Saca tu cuenta en GitHub, e instala el cliente en tu computadora a la vez.
RStudio te permitirá publicar tu trabajo final en un formato académico moderno usando markdown. Ten este tutorial a la mano todo el tiempo:https://rmarkdown.rstudio.com/index.html. Esto será tu último entregable.
Para mayores detalles, revisa el SILABO presente en Paideia y el campus virtual PUCP.
En esta sesión deseo que entiendas por qué es necesario ir más allá de la correlación (Pearson, Spearman, etc.) o las diferencias de valores centrales (t test, kruska wallis, etc.). Esta necesidad de ir más allá trae una técnica conocida como la regresión.
La sesión va en dos partes:
Trabajemos con estos datos:
library(rio)
topCloud='https://github.com/PoliticayGobiernoPUCP/'
courseCloud='estadistica_anapol2/raw/master/DATA/'
fileCloud='hsb_ok.xlsx'
hsb=import(paste0(topCloud,courseCloud,fileCloud))
Antes de correr cualquier análisis estadístico, debes revisar como el tipo de datos que tu archivo trae es reconocido por el R:
str(hsb)
## 'data.frame': 600 obs. of 15 variables:
## $ ID : chr "1" "2" "3" "4" ...
## $ SEX : num 2 1 2 2 2 1 1 2 1 2 ...
## $ RACE : num 2 2 2 2 2 2 2 2 2 2 ...
## $ SES : num 1 1 1 2 2 2 1 1 2 1 ...
## $ SCTYP : num 1 1 1 1 1 1 1 1 1 1 ...
## $ HSP : num 3 2 2 3 3 2 1 1 1 1 ...
## $ LOCUS : num 0.29 -0.42 0.71 0.06 0.22 0.46 0.44 0.68 0.06 0.05 ...
## $ CONCPT: num 0.88 0.03 0.03 0.03 -0.28 0.03 -0.47 0.25 0.56 0.15 ...
## $ MOT : num 0.67 0.33 0.67 0 0 0 0.33 1 0.33 1 ...
## $ CAR : num 10 2 9 15 1 11 10 9 9 11 ...
## $ RDG : num 33.6 46.9 41.6 38.9 36.3 49.5 62.7 44.2 46.9 44.2 ...
## $ WRTG : num 43.7 35.9 59.3 41.1 48.9 46.3 64.5 51.5 41.1 49.5 ...
## $ MATH : num 40.2 41.9 41.9 32.7 39.5 46.2 48 36.9 45.3 40.5 ...
## $ SCI : num 39 36.3 44.4 41.7 41.7 41.7 63.4 49.8 47.1 39 ...
## $ CIV : num 40.6 45.6 45.6 40.6 45.6 35.6 55.6 55.6 55.6 50.6 ...
Casi todo salió numérico, pero no sabremos qué ajustar si no leemos el manual metodológico o el diccionario de datos o la metadata disponible.
De ahi que debemos pre procesar:
categoricals=c("SEX","RACE","SES","SCTYP","HSP","CAR")
hsb[,categoricals]=lapply(hsb[,categoricals], as.factor)
# nominales
hsb$SEX=factor(hsb$SEX,
levels=c(1,2),
labels=c("Male","Female"))
hsb$RACE=factor(hsb$RACE,
levels=c(1,2,3,4),
labels=c("Hispanic","Asian","Black","White"))
hsb$HSP=factor(hsb$HSP,
levels=c(1,2,3),labels=c("General","Academic","Vocational"))
hsb$SCTYP=factor(hsb$SCTYP,
levels=c(1,2),
labels=c("Public","Private"))
# a ordinal:
hsb$SES=ordered(hsb$SES,
levels=c(1,2,3),
labels=c("Low","Medium","High" ))
Asumamos que nuestra variable de interés es el desempeño en matemáticas; así, nuestra variable dependiente está representada por la variable MATH.
A partir de ahí, consideremos que nos interesa saber la posible relación que pueda tener la variable que ha medido el desempeño en escritura; así, una variable independiente sería la representada por la variable WRTG. Hasta ahora sabemos que como son dos variables de tipo numérico debemos usar una correlación. La gráfica de correlación es esta:
library(ggplot2)
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
scatter = base + geom_point()
scatter
Vemos que hay una aparente relación. Calculemos los indices de correlación:
f1=formula(~MATH + WRTG)
# camino parametrico
pearsonf1=cor.test(f1,data=hsb)[c('estimate','p.value')]
# camino no parametrico
spearmanf1=cor.test(f1,data=hsb,method='spearman')[c('estimate','p.value')]
Asumiendo un camino paramétrico, podemos pedir el coeficiente de Pearson, el cuál al ser calculado obtenemos 0.6326664 (con p-value= 0).
Si hubieramos seguido una ruta no paramétrica, informaríamos el coeficiente de Spearman:0.6415126 (con p-value= 0).
Ahora, consideremos que nos interesa saber a la vez la posible relación que pueda tener la variable que ha medido el desempeño en ciencias; así otra variable independiente sería la representada por la variable SCI. Como es otra variable numérica no podemos calcular la correlación de tres variables, pero podemos tratar de verlo visualmente:
library(scatterplot3d)
scatterplot3d(hsb[,c('SCI','WRTG','MATH')])
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SCI))
Calculemos las correlaciones:
f2=formula(~MATH+SCI)
# camino parametrico
pearsonf2=cor.test(f2,data=hsb)[c('estimate','p.value')]
# camino no parametrico
spearmanf2=cor.test(f2,data=hsb, method='spearman')[c('estimate','p.value')]
Podríamos calcular la correlación de SCI con MATH, obteniendo el Pearson (0.6495261, p-value= 0) y el Spearman (0.6551515,p-value= 0).
Visualmente vemos relación, pero no tenemos un coeficiente para medir ello.
Finalmente, ¿si quisiéramos ver si el sexo tiene algun rol en todo esto? Como ésta es una variable categórica y dicotómica, lo primero que puede venir a la mente es esta gráfica:
base=ggplot(data=hsb, aes(x=as.factor(SEX), y=MATH))
base + geom_boxplot(notch = T) + geom_jitter(color="black", size=0.4, alpha=0.9)
Los boxplots tienen un notch flanqueando a la mediana, para sugerir igualdad de medianas si éstos se intersectan; de ahi que parece no haber diferencia sustantiva entre hombres y mujeres en cuanto a su desempeño en MATH.
Una alternativa al boxplot seria las barras de error:
library(ggpubr)
ggerrorplot(data=hsb, x = "SEX", y = "MATH")
En este último caso, si las lineas (denotado por las barras de error de la media) se intersectan, sugeriria que los valores medios (en este caso la media) podrian ser iguales.
Verificar si hay o no igualdad entre distribuciones depende si las variables se distribuyen o no de manera normal:
library(ggplot2)
ggplot(hsb,aes(x=MATH)) +
geom_histogram(aes(y = ..density..),bins = 20, fill='green') +
stat_function(fun = dnorm, colour = "red",
args = list(mean = mean(hsb$MATH, na.rm = TRUE),
sd = sd(hsb$MATH, na.rm = TRUE))) +
facet_grid(~SEX) +
coord_flip()
Nota que los histogramas de la data real tienen encima la curva normal que idealmente tendría esa data. La lejanía entre ellos, sugeriría no normalidad.
Se suele usar un qqplot para explorar la presencia/ausencia de normalidad:
# se sugiere normalidad si los puntos no se alejan de la diagonal.
library(ggpubr)
ggqqplot(data=hsb,x="MATH") + facet_grid(. ~ SEX)
Como ello no es fácil de discernir visualmente, tenemos por costumbre calcular algun coeficiente, como el Shapiro-Wilk:
library(knitr)
library(magrittr)
library(kableExtra)
f4=formula(MATH~SEX)
tablag= aggregate(f4, hsb,
FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})
# para que se vea mejor:
shapiroTest=as.data.frame(tablag[,2])
names(shapiroTest)=c("W","Prob")
kable(cbind(tablag[1],shapiroTest))%>%
kable_styling(bootstrap_options = c("striped", "hover"), full_width = F, position = "left")
| SEX | W | Prob |
|---|---|---|
| Male | 0.9837903 | 0.0034565 |
| Female | 0.9790040 | 0.0001031 |
Esto nos sugiere un camino no paramétrico para ver la diferencia de valores medios, por lo que deberiamos usar la prueba de Mann-Whitney en vez de la prueba t para testaer la relación entre ambas.
tf4=t.test(f4,data=hsb)[c('estimate','p.value')]
wilcoxf4=wilcox.test(f4,data=hsb)['p.value']
La prueba no paramétrica no rechazaría la igualdad de valores medios (Mann-Whitney con p valor = 0.3085543).
base=ggplot(data=hsb, aes(x=WRTG, y=MATH))
base + geom_point(aes(color = SEX))
base + geom_point(aes(size = SCI, color=SEX))
Otra alternativa puede ser:
base + geom_point(aes(color = SCI)) + facet_grid(~SEX)
Y claro:
paleta <- c("coral1","cyan" )
colors <- paleta[as.numeric(hsb$SEX)]
scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)
En todos los gráficos vemos que los hombres y las mujeres están distribuidos por todo el gráfico, lo cual nos sugiere que no hay diferencias aun en dimensiones mayores a dos. Sin embargo, no tenemos una medida de cuanto cada uno afecta a nuestra dependiente.
De ahi que necesitamos la regresión.
La regresión es una técnica donde hay que definir una variable dependiente y una o más independientes. Las independientes pueden tener rol predictor, dependiendo del diseño de investigación, aunque por defecto tiene un rol explicativo.
La regresión sí quiere informar cuánto una variable (independiente) puede explicar la variación de otra (dependiente), de ahí que es una técnica para probar hipótesis direccionales o asimétricas (las correlaciones tiene hipótesis simétricas).
La regresión busca proponer un modelo, es decir una ecuación, que recoja como una variable explicaría a otra.Para nuestro caso la variable dependiente es MATH, y tendriamos hasta ahora tres modelos:
modelo1=formula(MATH~WRTG)
modelo2=formula(MATH ~ WRTG + SCI)
modelo3= formula(MATH ~ WRTG + SCI + SEX)
Por ejemplo, para la hipótesis ‘el nivel de desempeño en escritura afecta el desempeño en matemáticas’, la regresión arrojaría este resultado:
O en mejor version con ayuda de stargazer:
library(stargazer)
reg1=lm(modelo1,data=hsb)
stargazer(reg1,type = "html",intercept.bottom = FALSE)
| Dependent variable: | |
| MATH | |
| Constant | 19.769*** |
| (1.633) | |
| WRTG | 0.612*** |
| (0.031) | |
| Observations | 600 |
| R2 | 0.400 |
| Adjusted R2 | 0.399 |
| Residual Std. Error | 7.297 (df = 598) |
| F Statistic | 399.110*** (df = 1; 598) |
| Note: | p<0.1; p<0.05; p<0.01 |
Aquí ya sabemos algo interesante, primero que WRTG tiene efecto, pues es significativo (indicado por los dos asteriscos); segundo, que ese efecto es directo, pues el coeficiente calculado es positivo; y tercero que la magnitud de ese efecto es 0.612, lo que indica cuanto aumenta MATH en promedio cuando WRTG se incremente en una unidad.
Esto es información suficiente para representar esa relación con una ecuación. Como la ecuación sólo tiene una variable independiente podemos producir una recta sobre el gráfico de correlación:
ggplot(hsb, aes(x=WRTG, y=MATH)) +
geom_point()+
geom_smooth(method=lm)
Esa recta podemos representarla así:
\[ MATH= 19.7690323 + 0.6123904 \cdot WRTG + \epsilon\]
El Y verdadero es MATH, pero la regresión produce un \(\hat{MATH}\) estimado, de ahi la presencia del \(\epsilon\). Justamente el R cuadrado ajustado (0.4002668) nos brinda un porcentaje (multiplicalo por 100) que da una pista de nuestra cercanía a una situación perfecta (cuando vale 1).
Y sí queremos ver el efecto de SCI?
reg2=lm(modelo2,data=hsb)
stargazer(reg2,type = "html",intercept.bottom = FALSE)
| Dependent variable: | |
| MATH | |
| Constant | 10.629*** |
| (1.630) | |
| WRTG | 0.377*** |
| (0.033) | |
| SCI | 0.415*** |
| (0.033) | |
| Observations | 600 |
| R2 | 0.524 |
| Adjusted R2 | 0.523 |
| Residual Std. Error | 6.505 (df = 597) |
| F Statistic | 328.846*** (df = 2; 597) |
| Note: | p<0.1; p<0.05; p<0.01 |
En este caso, la regresión tendrá una formula con dos variables explicando la dependiente, así que en vez de producir una línea buscará producir un plano:
library(scatterplot3d)
G <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')])
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)
Este plano podemos representarlo así:
\[ MATH= 10.6285904 + 0.3765304 \cdot WRTG + 0.4152733 \cdot SCI + \epsilon\]
Nuevamente, el Y verdadero es MATH, pero la regresión produce un \(\hat{MATH}\) estimado en forma de plano. De igual manera el R cuadrado ajustado (0.5241861) nos da una pista de nuestra lejanía a una situación perfecta.
Es clave darse cuenta de otro detalle, que el coeficiente de WRTG ha variado en la fórmula ahora que está presente SCI ¿Por qué sucede esto? Veamoslo así: en el primer caso, WRTG y \(\epsilon\) buscaban representar la variabilidad en MATH, y ahora, en el segundo caso, viene SCI para mejorar esa explicación; es así que el peso de la explicación ahora se recalcula y el coeficiente de WRTG deja de explicar lo que le corresponde a SCI, y \(\epsilon\) también le entrega algo a SCI.
Como \(\epsilon\) no tiene coeficiente, representamos su variación usando el error típico de los residuos o residual standard error (RSE). Nótese que éste ha variado de un modelo ha otro, ahora tenemos un RSE menor. Aquí vale la pena preguntarse si esta disminución del error es significativa, obteniendo:
tanova=anova(reg1,reg2)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")
| Res.Df | RSS | Df | Sum of Sq | F | Pr(> F) | |
| 1 | 598 | 31,842.070 | ||||
| 2 | 597 | 25,262.730 | 1 | 6,579.335 | 155.481 | 0 |
La comparación de modelos usando la tabla de análisis de varianza (anova) propone como hipótesis nula que los modelos no difieren (no se ha reducido el error al pasar de un modelo a otro). Como la comparación es significativa (vea el P), rechazamos igualdad de modelos: el modelo 2 sí reduce el error al incluir una variable más.
Finalmente, veamos el rol de sexo:
reg3=lm(modelo3,data=hsb)
stargazer(reg3,type = "html",intercept.bottom = FALSE)
| Dependent variable: | |
| MATH | |
| Constant | 11.306*** |
| (1.630) | |
| WRTG | 0.424*** |
| (0.036) | |
| SCI | 0.375*** |
| (0.035) | |
| SEXFemale | -1.922*** |
| (0.582) | |
| Observations | 600 |
| R2 | 0.533 |
| Adjusted R2 | 0.530 |
| Residual Std. Error | 6.452 (df = 596) |
| F Statistic | 226.510*** (df = 3; 596) |
| Note: | p<0.1; p<0.05; p<0.01 |
Aunque no podemos graficar cuatro coordenadas, podemos usar elementos visuales:
library(scatterplot3d)
colors <- paleta[as.numeric(hsb$SEX)]
G <- scatterplot3d(hsb[,c('SCI','WRTG','MATH')],color=colors)
G$plane3d(reg2, draw_polygon = TRUE, draw_lines = FALSE)
Nuestra nueva ecuación sería:
\[ MATH= 11.306349 + 0.4235786 \cdot WRTG + 0.3748025 \cdot SCI + -1.9219692 \cdot SEX + \epsilon\]
Nuevamente podemos ver si añadir SEXO en este modelo representa una mejora al anterior:
tanova=anova(reg1,reg2,reg3)
stargazer(tanova,type = 'html',summary = F,title = "Table de Análisis de Varianza")
| Res.Df | RSS | Df | Sum of Sq | F | Pr(> F) | |
| 1 | 598 | 31,842.070 | ||||
| 2 | 597 | 25,262.730 | 1 | 6,579.335 | 158.063 | 0 |
| 3 | 596 | 24,808.420 | 1 | 454.308 | 10.914 | 0.001 |
Finalmente, podemos resumir todo en esta tabla:
library(stargazer)
stargazer(reg1,reg2,reg3, type = "html", title = "Modelos planteadas",digits = 2, single.row = F,no.space = F,intercept.bottom = FALSE,
dep.var.caption="Variable dependiente:",
dep.var.labels="Desempeño en Matemáticas",
covariate.labels=c("Constante","Desempeño en Escritura","Desempeño en Ciencias","SEXO (mujer)"),
keep.stat = c("n","adj.rsq","ser"),df = F,
notes.label = "Notas:")
| Variable dependiente: | |||
| Desempeño en Matemáticas | |||
| (1) | (2) | (3) | |
| Constante | 19.77*** | 10.63*** | 11.31*** |
| (1.63) | (1.63) | (1.63) | |
| Desempeño en Escritura | 0.61*** | 0.38*** | 0.42*** |
| (0.03) | (0.03) | (0.04) | |
| Desempeño en Ciencias | 0.42*** | 0.37*** | |
| (0.03) | (0.04) | ||
| SEXO (mujer) | -1.92*** | ||
| (0.58) | |||
| Observations | 600 | 600 | 600 |
| Adjusted R2 | 0.40 | 0.52 | 0.53 |
| Residual Std. Error | 7.30 | 6.51 | 6.45 |
| Notas: | p<0.1; p<0.05; p<0.01 | ||
Gráficamente:
library(ggplot2)
library(sjPlot)
## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
plot_models(reg1,reg2,reg3,vline.color = "grey",m.labels=c("Modelo 1","Modelo 2","Modelo 3"))
Como vemos, el propósito del último gráfico es mostrar que dado que ninguna de las variables (y sus intervalos de confianza) tocan el valor cero (que significa que la variable no tiene efecto en la dependiente).